Zimbabwe Soybean Trials
Loading data
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsdata <- read.csv("data/Zimbabwe.csv",header = T,sep = ";")
# colnames(data)
# head(data)
pheno<- data[, c(4, 7, 16,21,25,27,38,41,46,63,64,65,66,68)]
table(pheno$LOCATION_NAME)##
## Banket_SCL Bindura_SCL Chisumbanje_CBI Goromonzi_SCL Gwebi_CBI
## 288 288 270 558 198
## Harare_CBI Harare_KKR Harare_SCL Kadoma_CBI Mazowe_TCK
## 468 90 288 180 36
## Mthampden_TCK Mutare_SCL Shamva_CBI Stapleford_SCL
## 90 378 360 468## [1] "CROP_SEASON" "TRIAL_NAME" "LOCATION_NAME" "REP_NO"
## [5] "Crop_Variety" "Control" "PL_Lodging" "R8_PL_Height"
## [9] "GRAIN_YLD_13PCT" "PLANT_DATE" "HARV_DATE" "SITE_LAT"
## [13] "SITE_LONG" "IRRIGATION"pheno$geno <- gsub(" ", "_",pheno$Crop_Variety)
pheno$geno <- factor(pheno$geno)
# unique(pheno$IRRIGATION)
levels(pheno$geno) <- paste("G", 1:length(levels(pheno$geno)), sep = "_")
# table(pheno$Crop_Variety)
# table(pheno$geno)
pheno[which(pheno$HARV_DATE==""),"LOCATION_NAME"]## [1] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [5] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [9] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [13] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [17] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [21] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [25] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [29] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [33] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [37] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [41] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [45] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [49] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [53] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [57] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [61] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [65] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [69] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [73] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [77] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [81] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [85] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [89] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [93] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [97] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [101] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [105] "Banket_SCL" "Banket_SCL" "Banket_SCL" "Banket_SCL"
## [109] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [113] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [117] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [121] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [125] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [129] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [133] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [137] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [141] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [145] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [149] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [153] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [157] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [161] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [165] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [169] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [173] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [177] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [181] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [185] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [189] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [193] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [197] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [201] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [205] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [209] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [213] "Bindura_SCL" "Bindura_SCL" "Bindura_SCL" "Bindura_SCL"
## [217] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [221] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [225] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [229] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [233] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [237] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [241] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [245] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [249] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [253] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [257] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [261] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [265] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [269] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [273] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [277] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [281] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [285] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [289] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [293] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [297] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [301] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [305] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [309] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [313] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [317] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [321] "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL" "Goromonzi_SCL"
## [325] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [329] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [333] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [337] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [341] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [345] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [349] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [353] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [357] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [361] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [365] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [369] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [373] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [377] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [381] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [385] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [389] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [393] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [397] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [401] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [405] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [409] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [413] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [417] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [421] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [425] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [429] "Harare_SCL" "Harare_SCL" "Harare_SCL" "Harare_SCL"
## [433] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [437] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [441] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [445] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [449] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [453] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [457] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [461] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [465] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [469] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [473] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [477] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [481] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [485] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [489] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [493] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [497] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [501] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [505] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [509] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [513] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [517] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [521] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [525] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [529] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [533] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [537] "Mutare_SCL" "Mutare_SCL" "Mutare_SCL" "Mutare_SCL"
## [541] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [545] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [549] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [553] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [557] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [561] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [565] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [569] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [573] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [577] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [581] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [585] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [589] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [593] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [597] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [601] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [605] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [609] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [613] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [617] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [621] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [625] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [629] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [633] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [637] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [641] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [645] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [649] "Harare_CBI"pheno$DTH <- as.Date(as.character(pheno$HARV_DATE), format="%m/%d/%Y")-
as.Date(as.character(pheno$PLANT_DATE), format="%m/%d/%Y")-5
# dim(pheno)
# colnames(pheno)
# head(pheno)
# str(pheno)
pheno$GRAIN_YLD_13PCT <- gsub(",", ".", pheno$GRAIN_YLD_13PCT)
trial <- gsub("_.*", "", pheno$LOCATION_NAME)
pheno$trial <- paste0(trial,"_", pheno$CROP_SEASON)
unique(pheno$TRIAL_NAME) # 41 trials## [1] "S18_2_ZW_SCL5" "S18_2_ZW_SCL11" "S18_2_ZW_CBI4" "S18_2_ZW_SCL7"
## [5] "S18_2_ZW_CBI1" "S18_2_ZW_SCL6" "S18_2_ZW_CBI2" "S18_2_ZW_SCL10"
## [9] "S18_2_ZW_CBI9" "S18_2_ZW_SCL8" "S19_1_ZW_CBI2" "S19_1_ZW_CBI4"
## [13] "S19_2_ZW_SCL5" "S19_2_ZW_CBI4" "S19_2_ZW_SCL1" "S19_2_ZW_CBI1"
## [17] "S19_2_ZW_SCL2" "S19_2_ZW_SCL6" "S19_2_ZW_CBI5" "S19_2_ZW_SCL3"
## [21] "S20_2_ZW_SCL10" "S20_2_ZW_SCL11" "S20_2_ZW_SCL7" "S20_2_ZW_CBI2"
## [25] "S20_2_ZW_CBI1" "S20_2_ZW_SCL8" "S20_2_ZW_SCL12" "S20_2_ZW_SCL9"
## [29] "S21_2_ZW_SCL3" "S21_2_ZW_SCL2" "S21_2_ZW_CBI2" "S21_2_ZW_CBI1"
## [33] "S21_2_ZW_KKR1" "S21_2_ZW_CBI3" "S21_2_ZW_SCL1" "S22_2_ZW_SCL3"
## [37] "S22_2_ZW_CBI1" "S22_2_ZW_TCK1" "S22_2_ZW_CBI3" "S22_2_ZW_SCL4"
## [41] "S23_2_ZW_SCL1" "S23_2_ZW_TCK1" "S23_2_ZW_SCL2"## [1] "Banket_SCL" "Bindura_SCL" "Chisumbanje_CBI" "Goromonzi_SCL"
## [5] "Harare_CBI" "Harare_SCL" "Kadoma_CBI" "Mutare_SCL"
## [9] "Shamva_CBI" "Stapleford_SCL" "Gwebi_CBI" "Harare_KKR"
## [13] "Mazowe_TCK" "Mthampden_TCK"## [1] 2018 2019 2020 2021 2022 2023## [1] 98## CROP_SEASON TRIAL_NAME LOCATION_NAME REP_NO Crop_Variety Control
## 1 2018 S18_2_ZW_SCL5 Banket_SCL 1 S1079/6/7 No
## 2 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 1987-62F No
## 3 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 1991-22F No
## 4 2018 S18_2_ZW_SCL5 Banket_SCL 1 Lukanga No
## 5 2018 S18_2_ZW_SCL5 Banket_SCL 1 S1146/5/25 No
## 6 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 2014-16FM No
## PL_Lodging R8_PL_Height GRAIN_YLD_13PCT PLANT_DATE HARV_DATE SITE_LAT
## 1 1 105 4680.55 12/04/2018 4/26/2019 -176,425
## 2 5 95 2236.97 12/04/2018 4/26/2019 -176,425
## 3 5 90 1866.21 12/04/2018 4/26/2019 -176,425
## 4 1 100 3723.68 12/04/2018 4/26/2019 -176,425
## 5 1 90 3185.52 12/04/2018 4/26/2019 -176,425
## 6 5 100 2851.49 12/04/2018 4/26/2019 -176,425
## SITE_LONG IRRIGATION geno DTH trial
## 1 30,399,997 Rainfed G_25 138 days Banket_2018
## 2 30,399,997 Rainfed G_56 138 days Banket_2018
## 3 30,399,997 Rainfed G_57 138 days Banket_2018
## 4 30,399,997 Rainfed G_8 138 days Banket_2018
## 5 30,399,997 Rainfed G_27 138 days Banket_2018
## 6 30,399,997 Rainfed G_69 138 days Banket_2018pheno <- pheno |> rename("Plant Height"=8, "Grain Yield"=9, "Plant lodging"=7, "rep"=4)
traits <- c(colnames(pheno)[c(7,8,9)])
traits## [1] "Plant lodging" "Plant Height" "Grain Yield"## Warning in lapply(pheno[, traits], as.numeric): NAs introduced by coercion## [1] "CROP_SEASON" "TRIAL_NAME" "LOCATION_NAME" "rep"
## [5] "Crop_Variety" "Control" "IRRIGATION" "geno"
## [9] "DTH" "trial"## [1] "Plant lodging" "Plant Height" "Grain Yield"## CROP_SEASON TRIAL_NAME LOCATION_NAME rep Crop_Variety Control
## 1 2018 S18_2_ZW_SCL5 Banket_SCL 1 S1079/6/7 No
## 2 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 1987-62F No
## 3 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 1991-22F No
## 4 2018 S18_2_ZW_SCL5 Banket_SCL 1 Lukanga No
## 5 2018 S18_2_ZW_SCL5 Banket_SCL 1 S1146/5/25 No
## 6 2018 S18_2_ZW_SCL5 Banket_SCL 1 TGx 2014-16FM No
## Plant lodging Plant Height Grain Yield PLANT_DATE HARV_DATE SITE_LAT
## 1 1 105 4680.55 12/04/2018 4/26/2019 -176,425
## 2 5 95 2236.97 12/04/2018 4/26/2019 -176,425
## 3 5 90 1866.21 12/04/2018 4/26/2019 -176,425
## 4 1 100 3723.68 12/04/2018 4/26/2019 -176,425
## 5 1 90 3185.52 12/04/2018 4/26/2019 -176,425
## 6 5 100 2851.49 12/04/2018 4/26/2019 -176,425
## SITE_LONG IRRIGATION geno DTH trial
## 1 30,399,997 Rainfed G_25 138 Banket_2018
## 2 30,399,997 Rainfed G_56 138 Banket_2018
## 3 30,399,997 Rainfed G_57 138 Banket_2018
## 4 30,399,997 Rainfed G_8 138 Banket_2018
## 5 30,399,997 Rainfed G_27 138 Banket_2018
## 6 30,399,997 Rainfed G_69 138 Banket_2018#Irrigation system
# table(pheno$trial,pheno$IRRIGATION, useNA = "ifany")
# unique(pheno$IRRIGATION)
pheno[which(pheno$trial=="Banket_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Banket_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Bindura_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Chisumbanje_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Goromonzi_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Goromonzi_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Gwebi_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Mutare_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial==" Mutare_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Stapleford_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Stapleford_2020"),"IRRIGATION"] <- "Rainfed"
# table(pheno$IRRIGATION)
#
#
# levels(pheno$IRRIGATION)Data visualization
pheno_long <- reshape2::melt(pheno, measure.vars = traits, variable.name = "trait") #transformando o pheno em formato longo
# head(pheno_long)
ggplot(pheno_long, aes(x = trait, y = value, color = trait)) + geom_boxplot(outlier.shape = NA)+
geom_point(alpha = 0.2, position = "jitter", size = 1) + facet_grid(trait~LOCATION_NAME,
scales = "free_y") + theme_bw() + theme(axis.text.x = element_text(angle=90))## Warning: Removed 191 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 191 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Loading required package: ggpp## Registered S3 methods overwritten by 'ggpp':
## method from
## heightDetails.titleGrob ggplot2
## widthDetails.titleGrob ggplot2##
## Attaching package: 'ggpp'## The following object is masked from 'package:ggh4x':
##
## stat_centroid## The following object is masked from 'package:ggplot2':
##
## annotatedata_vis=ggplot(data = pheno_long, aes(x=trait, y= value, color=trait)) +
geom_boxplot() +
geom_point( position = "jitter", alpha = 0.5) + facet_nested(trait~LOCATION_NAME+CROP_SEASON, scales = "free_y") +
theme_bw();data_vis## Warning: Removed 191 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Removed 191 rows containing missing values or values outside the scale range
## (`geom_point()`).
locations <- pheno[,c("TRIAL_NAME","LOCATION_NAME","SITE_LAT","SITE_LONG")]
locations_unq <- locations |> distinct(SITE_LAT,.keep_all = TRUE)
locations_unq <- locations_unq[-21,]
smart_clean_coords <- function(x, decimal_places = 5) {
# Replace commas with dots (fix decimals)
x <- gsub(",", ".", x)
# Remove thousand separators (dots between digits)
x <- gsub("(?<=\\d)\\.(?=\\d)", "", x, perl = TRUE)
# Convert to numeric (handles scientific notation like 0.00E+00)
num <- as.numeric(x)
# Function to scale down numbers that are too large
scale_if_needed <- function(val) {
if (is.na(val)) return(NA)
# Latitude should be between -90 and 90
if (abs(val) > 90) {
# Keep dividing by 10 until in valid range
while (abs(val) > 90) {
val <- val / 10
}
return(val)
} else {
return(val) # Already valid, return as-is
}
}
# Apply scaling only where needed
num <- sapply(num, scale_if_needed)
# Round to specified decimal places
round(num, decimal_places)
}
locations_unq$SITE_LAT <- smart_clean_coords(locations_unq$SITE_LAT)
locations_unq$SITE_LONG <- smart_clean_coords(locations_unq$SITE_LONG)
write.csv(locations_unq,file = "data/loc_unq.csv",row.names = F)Mapping the trials
#Com ggplot2
# elevação
library(raster)
library(RColorBrewer)
library(ggplot2)
library(sf)
africa.sf<-st_read("data/afr_g2014_2013_0.shp")
# uf.sf<- brasil.sf[brasil.sf$UF %in% c('MT','GO','RS','PE','BA','ES','DF',
# 'RJ','AL','SE','MG','MS','SP','PR','SC'), ]
library(stringi)
africa.sf$ADM0_NAME <- stri_replace_all_fixed(
africa.sf$ADM0_NAME,
"C\xf4te d'Ivoire", # Literal special character
"Cote d'Ivoire"
)
# elevação
#Selecionando os estados
africa<- shapefile('data/afr_g2014_2013_0.shp')
africa
#brasil[LINHA , COLUNA]
uf<- africa[africa@data$ADM0_NAME%in% c("Zimbabwe"),] #Mesma indexacao de data.frame
# # uf
#
plot(uf, col='orange')
# plot(uf, col=c("red","blue"))
# text(uf,'UF', cex=1, halo=T)
#Cortando um raster para o estado
agro<- raster('data/agro.tif') # Importando o Raster de elevação
hwsd <- raster('data/hwsd.tif')
plot(agro)
plot(hwsd)
#for( i in 1:5){
hwsd_africa<- crop(africa,uf)
plot(hwsd_africa)
hwsd_africa<-mask(africa,uf)
plot(hwsd_africa)
palete.spec<- brewer.pal("Spectral",n=11)
#plot(elev_uf, col=rev(palete.spec))
writeRaster(hwsd,filename= "data/hwsd_map", overwrite=TRUE)
writeRaster(agro,filename= "data/agro_map", overwrite=TRUE)
WGS84<-CRS("+proj=longlat +datum=WGS84 +no_defs ") #
#pontos
pontos<- read.csv("data/loc_unq.csv", header=TRUE, sep=',') # se estiver usando outro OS
str(pontos)
# pontos[which(pontos$latitude < -23),]
# pontos[which(pontos$Instituicao=="CATI/Avare"),"regiao"]="SE"
# pontos[which(pontos$Instituicao=="DETEC/Itai"),"regiao"]="SE"
# pontos[which(pontos$cidade=="Campos Novos"),"regiao"]="SU"
# Apply to your data (5 decimal rounding)
# pontos$SITE_LAT <- smart_clean_coords(pontos$SITE_LAT)
# pontos$SITE_LONG <- smart_clean_coords(pontos$SITE_LONG)
# Check results
head(pontos)
tail(pontos)
coordinates(pontos)<- ~ SITE_LONG + SITE_LAT
str(pontos)
#plot(pontos,pch=20,cex=0.5, add=T)
crs(pontos)<-WGS84
# fv <- unique(pontos$regiao)
# (m <- match(pontos$regiao, fv))
### Plotando o mapa do Brasil
#
ggplot() +
geom_sf(data = africa.sf)
#scale_fill_discrete(type = c("blue",'yellow','red','magenta','green'))
### Plotando o mapa do Brasil e da regiao cortada
#
ggplot() +
geom_sf(data = africa.sf, fill="white")
pontos.sf
pontos
### Aproximando para a regi?o cortada
#
pontos@coords.nrs
(bbox_uf<-st_bbox(uf)) #Pegando os limites da regi?o cortada
ggplot() +
geom_sf(data = africa.sf, fill="white") +
#Código para definir os limites do mapa cood_sf()
coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"]))+
theme_light()
#
### Adicionando os rótulos de estado
#
ggplot() +
geom_sf(data = africa.sf, fill="grey") +
theme_light()
### Plotando pontos
#
(pontos.sf<-data.frame(pontos)) #retornando nossos pontos para tabela
# pontos.sf <-pontos.sf[-21,]
pontos.sf
# ggplot() +
# geom_sf(data = brasil.sf, fill="grey") +
# geom_sf(data = uf.sf,fill="white")+
# coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"])) +
# geom_sf_label(data = brasil.sf, aes(label = UF),fill=NA,label.size = NA)+
# #Código para plotagem de pontos (geom_point)
# geom_point(data = pontos.sf, mapping = aes(x = longitude, y = latitude, color=regiao),size =1) +
# scale_color_manual(breaks = c("NE", "CO", "SE","SU"),values=rainbow(4))+
# theme_light()+ guides(color = guide_legend(title = "Region"))+labs(y="Latitude",x="Longitude")
problems <- africa.sf$ADM0_NAME[!stringi::stri_enc_isascii(africa.sf$ADM0_NAME)]
if(length(problems) > 0) print(problems)
### com rasters
#
# install.packages("tmap")
library(stars)
elev_uf.stars<- read_stars("data/hwsd_map.grd",RasterIO = list(nBufXSize = 600, nBufYSize = 600))
# elev_uf.stars<- read_stars("data/agro_map.grd",RasterIO = list(nBufXSize = 600, nBufYSize = 600))
names(elev_uf.stars)<-"Elevation"
pontos
st_bbox(africa)
(bbox_uf<-st_bbox(uf))
(pontos.sf<-data.frame(pontos))
mapa=ggplot() +
geom_sf(data = africa.sf, fill="white") +
#Código para plotagem do rater geom_stars()+
coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"]),crs = "WGS84") +theme_light()+
geom_point(data = pontos.sf, mapping = aes(x = SITE_LONG, y = SITE_LAT),size =1,show.legend = T)+
geom_sf_label(data = africa.sf, aes(label = ADM0_NAME),fill=NA,label.size = NA)+theme(legend.position = "none",axis.text.x =element_text(angle = 90))+labs(x="",y="");mapa
africamap=ggplot() +
geom_sf(data = africa.sf, fill="white")+
scale_fill_distiller(palette = "RdYlGn",na.value = NA,name = "Altitude")+
geom_stars(data = elev_uf.stars,alpha =0.8)+
coord_sf(xlim =c(-25.35875,63.50265) ,ylim=c(-46.9813, 37.56095),crs = "WGS84")+labs(y="Latitude",x="Longitude")+
theme_light()+theme(legend.position = "none")+
geom_rect(data=africa.sf,
xmin = 25.23703,
ymin =-22.41774,
xmax = 33.05630,
ymax =-15.60883 ,
fill = NA,
colour = "black",
linewidth = 0.8
);africamap
library(cowplot)
map=africamap |> ggdraw()+
draw_plot(
{
mapa +
coord_sf(
xlim = c(25.23703, 33.05630),
ylim = c(-22.41774, -15.60883),
expand = FALSE) +
theme(legend.position = "none")
},
# The distance along a (0,1) x-axis to draw the left edge of the plot
x = 0.13,
# The distance along a (0,1) y-axis to draw the bottom edge of the plot
y = 0.05,
# The width and height of the plot expressed as proportion of the entire ggdraw object
width = 0.40,
height = 0.40);
map
ggsave(plot=map,device = "pdf",filename ="Plots/soybean_zimbabwe.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=map,device = "png",filename ="Plots/soybean_zimbabwe.png",dpi = "retina",height = 6,width = 8)
# Single-environment
analyses
These analyses had a purpose: to test the significance of the genetic variance in each environment. We used the following model:
\[ \mathbf{y} = \mathbf{X_1b} + \mathbf{X_2che} + \mathbf{Z_1g} + \mathbf{\epsilon} \]
where \(\mathbf{y}\) is the vector of phenotypic observations, \(\mathbf{b}\) is the vector of fixed effects of replication, \(\mathbf{che}\) is the vector of fixed effects of checks, \(\mathbf{g}\) is the vector of random effects of genotypes and \(\mathbf{\epsilon}\) the vector of errors. \(\mathbf{X_1}\), \(\mathbf{X_2}\), \(\mathbf{Z_1}\) $ are incidence matrix of \(\mathbf{b}\), \(\mathbf{g}\) effects respectively.
At each environment, we filtered out trials that do not have significance for genotype effect. The Likelihood Ratio Test (\(LRT\)) was computed as follows:
\[ LRT= −2 \times (Log𝐿 - Log L_𝑅) \]
where \(L\) is the maximum point of residual likelihood function of the complete model and \(L_R\) is the same for the reduced model, that is, without the effect to be tested. The LRT value was compared with a tabulated value based on the chi-square table, with one degree of freedom and 0.95 probability.
We calculated the broad-sense heritability (\(h^2\)):
\[ h^2 = \sigma^2g / \sigma^2g + \sigma^2e \]
where \(\sigma^2g\) is the is the additive genetic variance,\(\sigma^2g\) is non additive genetic variance and \(\sigma_e^2\) is the error variance.
STMT
source("https://raw.githubusercontent.com/saulo-chaves/May_b_useful/main/update_asreml.R")
#source("GEBV/codes/outliers.R")
library(asreml)
data.list = split(pheno, f = pheno$trial)
traits
stmt = list()
j="Banket_2018"
#j="2021NAC_ST01"
i="Plant lodging"
unique(pheno$trial)
for (j in names(data.list)) {
x = data.list[[j]]
cat("====> Environment:", j, fill = TRUE)
# A.env = full2sparse(K = amat[which(rownames(amat) %in%
# c(x$geno)),
# which(colnames(amat) %in%
# c(x$geno))])
# attr(A.env, "INVERSE") = FALSE
# H.env.van = full2sparse(K = H.inv.van[which(rownames(H.inv.van) %in%
# c(x$geno)),
# which(colnames(H.inv.van) %in%
# c(x$geno))])
#
# attr(H.env.van, "INVERSE") = TRUE
# x = x |> mutate(nad=geno) |> mutate_at(vars(env, rep, geno, row, col, check,nad), as.factor)
summar.mat = matrix(nrow = length(traits), ncol = 4,
dimnames = list(traits,
c("geno","error","H2", "CV")))
vc.env = list()
Blue.env = list()
for (i in traits) {
cat("====> Trait:", i, fill = TRUE)
if(all(is.na(x[,i]))) next
model.r = asreml(fixed = x[,i] ~ rep,
random = ~ geno,maxit = 50,
data = x, na.action = na.method(x = "include", y = "include"))
# colocar nugget and spline
model.f = asreml(fixed = x[,i] ~ rep+geno,
maxit = 50,
residual = ~units,
data = x, na.action = na.method(x = "include", y = "include"))
if(!model.r$converge) next
if(any(na.exclude(model.r$vparameters.pc > 1))) model.r = update.asreml(model.r)
# varcomp = summary(modelhvan)$varcomp
varcomp = summary(model.r)$varcomp
varcomp.df = data.frame(
effect = c('genotypic', "residual"),
component = c(round(
varcomp[grep('geno', rownames(varcomp)),1],4),
round(varcomp[grep('!R', rownames(varcomp)),1],4))
)
# if(!modelad$converge) next
#
# if(any(na.exclude(modelad$vparameters.pc > 1))) modelad = up.mod(modelad)
#
if(!model.f$converge) next
if(any(na.exclude(model.f$vparameters.pc > 1))) model = up.mod(model.f)
lrtest = as.data.frame(lrt(asreml(fixed = x[,i] ~ rep,
residual = ~units, maxit = 50,
data = x),model.r))
lrtest=up.mod(lrtest)
lrtest$effect = c("genotypic")
varcomp= full_join(varcomp.df, lrtest[,-1], by = "effect")
varcomp$signi = ifelse(varcomp$`Pr(Chisq)` <= 0.06, "*", "ns")
varcomp$signi2 = ifelse(varcomp$`Pr(Chisq)` <= 0.06,
paste(round(varcomp$component, 4), "*"),
paste(round(varcomp$component, 4), "ns"))
rownames(varcomp)=varcomp$effect
pred = predict(model.r, classify = "geno", sed = TRUE)
pred$sed= pred$sed[-which(pred$pvals$status=="Aliased"),-which(pred$pvals$status=="Aliased")]
MVdelta = mean((pred$sed^2)[upper.tri(pred$sed^2,diag = F)])
H2 = 1-(MVdelta/(2*varcomp[varcomp$effect == "genotypic","component"]))
h2 = vpredict(model.r, h2~(V1)/(V1+V2))
#h2ad = vpredict(modelad, h2~(V3)/(V3+V4))
CV = sqrt(varcomp[grep("r", rownames(varcomp)),"component"])/
mean(x[,i], na.rm = TRUE)*100
summar.mat[i,] = c(varcomp$signi2[1], round(varcomp$component[2],4),
round(h2$Estimate,4), round(CV,4))
pred=predict(model.f, classify = "geno", vcov = T)
pred$vcov= pred$vcov[-which(pred$pvals$status=="Aliased"),-which(pred$pvals$status=="Aliased")]
pred$pvals= pred$pvals[-which(pred$pvals$status=="Aliased"),]
pred$pvals$w=diag(solve(pred$vcov))
pred=pred$pvals
vc.env[[i]] = varcomp.df
Blue.env[[i]] = pred
}
stmt[[j]] = list(Blues = Blue.env,summary = summar.mat,vc=vc.env)
rm(summar.mat,GEBV.env, varcomp, lrtest, model,
h2, CV, x, H.inv.env)
}
saveRDS(stmt, file = "Saves/stmt.rds")
#rm(data.list)
#hist(data[which(data$env=="UGN2021ABI_ST02"),"rytha"])
#table(data[which(data$env=="UGN2021ABI_ST02"),"rytha"],useNA = "ifany")
#table(data[which(data$env=="UGN2020NAS_ST04"),"vir2"], useNA = "ifany")Results
stmt=readRDS("Saves/stmt.rds")
stvc=lapply(stmt, function(x)as.data.frame(cbind(x$summary)))
#stvc.df=lapply(stmti, function(x)as.data.frame(cbind(x$vc)))
#stmt[[1]][["Blues"]][[1]][,-c(3,4)]
stvc=do.call(rbind,stvc)
stvc$trait=rownames(stvc)
stvc$env=sub("\\..*","",stvc$trait)
stvc$trait = sub(".*\\.","",stvc$trait)
stvc$geno_sig=sub("^\\S+\\s+","",stvc$geno)
stvc$geno=sub("\\s.*","",stvc$geno)
stvc$geno= as.numeric(stvc$geno)
stvc$H2=as.numeric(stvc$H2)
stvc$error=as.numeric(stvc$error)
stvc$CV=as.numeric(stvc$CV)
library(tidyverse)
head(stvc)## geno error H2 CV trait
## Banket_2018.Plant lodging 0.7079 1.1989 0.3713 46.7037 Plant lodging
## Banket_2018.Plant Height 484.7230 123.1807 0.7974 10.5836 Plant Height
## Banket_2018.Grain Yield 165411.8384 324929.2135 0.3373 22.4875 Grain Yield
## Banket_2019.Plant lodging NA NA NA NA Plant lodging
## Banket_2019.Plant Height NA NA NA NA Plant Height
## Banket_2019.Grain Yield 184876.1880 295453.6310 0.3849 17.2950 Grain Yield
## env geno_sig
## Banket_2018.Plant lodging Banket_2018 *
## Banket_2018.Plant Height Banket_2018 *
## Banket_2018.Grain Yield Banket_2018 *
## Banket_2019.Plant lodging Banket_2019 <NA>
## Banket_2019.Plant Height Banket_2019 <NA>
## Banket_2019.Grain Yield Banket_2019 *stvcl= stvc |> select(geno, error, trait, env, geno_sig) |> pivot_longer(geno:error)
unique(stvc$env)## [1] "Banket_2018" "Banket_2019" "Banket_2020" "Bindura_2018"
## [5] "Bindura_2020" "Bindura_2021" "Chisumbanje_2018" "Chisumbanje_2019"
## [9] "Goromonzi_2018" "Goromonzi_2019" "Goromonzi_2020" "Goromonzi_2021"
## [13] "Goromonzi_2022" "Goromonzi_2023" "Gwebi_2020" "Gwebi_2021"
## [17] "Harare_2018" "Harare_2019" "Harare_2020" "Harare_2021"
## [21] "Harare_2022" "Kadoma_2018" "Kadoma_2021" "Mazowe_2022"
## [25] "Mthampden_2023" "Mutare_2018" "Mutare_2019" "Mutare_2020"
## [29] "Mutare_2021" "Shamva_2018" "Shamva_2019" "Shamva_2022"
## [33] "Stapleford_2018" "Stapleford_2019" "Stapleford_2020" "Stapleford_2022"
## [37] "Stapleford_2023"# cols.label=c("A20NAC", "A20NAS", "B20NAS", "A21ABI" ,"B21ABI","A21NAC", "B21NAC", "A21NAS", "B21NAS")
# names(cols.label)=c(unique(stvc$env))
stvcl|> ggplot(aes(x=name,y=value))+ geom_col(data= subset(stvcl, !grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(fill="*"), color="black", fill="#a8ddb5") + geom_col(data= subset(stvcl, grepl("error",stvcl$name)),color="black",fill="lightblue") + facet_grid(trait~env,scales = "free_y") + geom_text(data=subset(stvcl, !grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(label = geno_sig), size=6 ) + theme(legend.position = "top", ) + geom_col(data= subset(stvcl, grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(fill="ns"), fill="gray", color="black") + theme(legend.position = "top", ) + labs(x="Effects", y="Genotypic variance component") + geom_text(data=subset(stvcl, grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(label = "ns", vjust = -1.0), size=3 ) ## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Removed 17 rows containing missing values or values outside the scale range
## (`geom_col()`).## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_text()`).
stvc_bef <- stvc |> select(H2, CV, trait, env, geno_sig) |> pivot_longer(H2:CV) |>
ggplot(aes(x=env,y=value))+geom_col(aes(fill = trait), position = position_dodge())+ facet_grid(name~trait, scales = "free_y") + theme(legend.position = "none", axis.text.x = element_text(angle = 90)) + labs(x= "Environments", y="Parameters")
stvc$CV[stvc$CV == 948.7212] <- NA
stvc$CV[stvc$CV == 130.8657] <- NA
stvc$CV[stvc$CV == 129.7948] <- NA
stvc_aft <- stvc |> select(H2, CV, trait, env, geno_sig) |> pivot_longer(H2:CV) |>
ggplot(aes(x=env,y=value))+geom_col(aes(fill = trait), position = position_dodge())+ facet_grid(name~trait, scales = "free_y") + theme(legend.position = "none", axis.text.x = element_text(angle = 90)) + labs(x= "Environments", y="Parameters")
data.list = split(pheno, f = pheno$trial)
z = data.list[["Stapleford_2018"]]
z = data.list[["Harare_2021"]]
table(z$`Plant lodging`)##
## 0 1 2 3 4 5
## 91 38 6 22 12 11##
## 0 1 2 3 4 5
## 108 39 32 19 5 13##
## 0 1 2 3 4 5
## 9 16 10 14 6 35## geno error H2 CV trait
## Banket_2018.Plant lodging 0.7079 1.1989 0.3713 NA Plant lodging
## Banket_2018.Plant Height 484.7230 123.1807 0.7974 NA Plant Height
## Banket_2018.Grain Yield 165411.8384 324929.2135 0.3373 NA Grain Yield
## Banket_2019.Plant lodging NA NA NA NA Plant lodging
## Banket_2019.Plant Height NA NA NA NA Plant Height
## Banket_2019.Grain Yield 184876.1880 295453.6310 0.3849 NA Grain Yield
## Banket_2020.Plant lodging NA NA NA NA Plant lodging
## Banket_2020.Plant Height 199.5743 40.6721 0.8307 NA Plant Height
## Banket_2020.Grain Yield 300108.5338 891371.4859 0.2519 NA Grain Yield
## Bindura_2018.Plant lodging 1.7364 0.4628 0.7895 NA Plant lodging
## Bindura_2018.Plant Height 623.1939 88.5958 0.8755 NA Plant Height
## Bindura_2018.Grain Yield 330922.1419 172315.2347 0.6576 NA Grain Yield
## Bindura_2020.Plant lodging NA NA NA NA Plant lodging
## Bindura_2020.Plant Height 33.9428 285.5884 0.1062 NA Plant Height
## Bindura_2020.Grain Yield 189852.6266 697750.7882 0.2139 NA Grain Yield
## Bindura_2021.Plant lodging NA NA NA NA Plant lodging
## Bindura_2021.Plant Height NA NA NA NA Plant Height
## Bindura_2021.Grain Yield 365074.5146 250284.6771 0.5933 NA Grain Yield
## Chisumbanje_2018.Plant lodging 1.5659 0.5192 0.7510 NA Plant lodging
## Chisumbanje_2018.Plant Height 845.3728 124.2484 0.8719 NA Plant Height
## Chisumbanje_2018.Grain Yield 259732.3391 230272.7664 0.5301 NA Grain Yield
## Chisumbanje_2019.Plant lodging 0.7430 0.9201 0.4468 NA Plant lodging
## Chisumbanje_2019.Plant Height 158.9274 118.1124 0.5737 NA Plant Height
## Chisumbanje_2019.Grain Yield 9787.6001 272120.6313 0.0347 NA Grain Yield
## Goromonzi_2018.Plant lodging NA NA NA NA Plant lodging
## Goromonzi_2018.Plant Height 539.0736 67.5482 0.8886 NA Plant Height
## Goromonzi_2018.Grain Yield 184716.1170 168639.1436 0.5227 NA Grain Yield
## Goromonzi_2019.Plant lodging 0.2935 0.2571 0.5331 NA Plant lodging
## Goromonzi_2019.Plant Height 283.4581 90.1553 0.7587 NA Plant Height
## Goromonzi_2019.Grain Yield 301347.8832 342784.4546 0.4678 NA Grain Yield
## Goromonzi_2020.Plant lodging NA NA NA NA Plant lodging
## Goromonzi_2020.Plant Height 140.9136 50.0625 0.7379 NA Plant Height
## Goromonzi_2020.Grain Yield 281956.8209 414415.4873 0.4049 NA Grain Yield
## Goromonzi_2021.Plant lodging NA NA NA NA Plant lodging
## Goromonzi_2021.Plant Height 193.6497 53.6208 0.7831 NA Plant Height
## Goromonzi_2021.Grain Yield 238728.0291 87261.9853 0.7323 NA Grain Yield
## Goromonzi_2022.Plant lodging 0.4584 0.8481 0.3508 NA Plant lodging
## Goromonzi_2022.Plant Height 338.2353 207.7936 0.6194 NA Plant Height
## Goromonzi_2022.Grain Yield 100029.8381 319937.0838 0.2382 NA Grain Yield
## Goromonzi_2023.Plant lodging NA NA NA NA Plant lodging
## Goromonzi_2023.Plant Height 22.6735 166.5719 0.1198 NA Plant Height
## Goromonzi_2023.Grain Yield 538135.4862 256713.1407 0.6770 NA Grain Yield
## Gwebi_2020.Plant lodging 1.1696 0.7089 0.6226 NA Plant lodging
## Gwebi_2020.Plant Height 97.8728 110.9833 0.4686 NA Plant Height
## Gwebi_2020.Grain Yield 11419.4391 471465.2680 0.0236 NA Grain Yield
## Gwebi_2021.Plant lodging 2.1307 0.2502 0.8949 NA Plant lodging
## Gwebi_2021.Plant Height 218.4241 37.0886 0.8548 NA Plant Height
## Gwebi_2021.Grain Yield 385098.9591 468916.0707 0.4509 NA Grain Yield
## Harare_2018.Plant lodging 0.7686 0.5908 0.5654 NA Plant lodging
## Harare_2018.Plant Height 527.0460 98.0460 0.8431 NA Plant Height
## Harare_2018.Grain Yield 106664.3020 326694.1673 0.2461 NA Grain Yield
## Harare_2019.Plant lodging 0.3409 1.4088 0.1948 NA Plant lodging
## Harare_2019.Plant Height 80.7317 228.7499 0.2609 NA Plant Height
## Harare_2019.Grain Yield 0.1686 1666412.8585 0.0000 NA Grain Yield
## Harare_2020.Plant lodging 0.0000 2.1674 0.0000 NA Plant lodging
## Harare_2020.Plant Height 99.2981 148.1606 0.4013 NA Plant Height
## Harare_2020.Grain Yield 754984.8867 810821.7839 0.4822 NA Grain Yield
## Harare_2021.Plant lodging 0.0586 2.5351 0.0226 NA Plant lodging
## Harare_2021.Plant Height 106.7914 30.6941 0.7767 NA Plant Height
## Harare_2021.Grain Yield 322033.7747 460802.9537 0.4114 NA Grain Yield
## Harare_2022.Plant lodging 1.3000 1.9318 0.4023 NA Plant lodging
## Harare_2022.Plant Height 16.6199 163.4875 0.0923 NA Plant Height
## Harare_2022.Grain Yield 42130.5542 471560.2576 0.0820 NA Grain Yield
## Kadoma_2018.Plant lodging 0.8203 1.1468 0.4170 NA Plant lodging
## Kadoma_2018.Plant Height 57.1665 149.5267 0.2766 NA Plant Height
## Kadoma_2018.Grain Yield 447777.9957 1003626.5721 0.3085 NA Grain Yield
## Kadoma_2021.Plant lodging 0.1188 1.9046 0.0587 NA Plant lodging
## Kadoma_2021.Plant Height 35.3694 311.7297 0.1019 NA Plant Height
## Kadoma_2021.Grain Yield 33215.5894 2089801.1610 0.0156 NA Grain Yield
## Mazowe_2022.Plant lodging 0.0152 0.4066 0.0359 NA Plant lodging
## Mazowe_2022.Plant Height 219.8637 35.4920 0.8610 NA Plant Height
## Mazowe_2022.Grain Yield 42707.5208 161563.0444 0.2091 NA Grain Yield
## Mthampden_2023.Plant lodging 0.2155 0.3911 0.3552 NA Plant lodging
## Mthampden_2023.Plant Height 94.3963 51.7369 0.6460 NA Plant Height
## Mthampden_2023.Grain Yield 25746.4389 157702.2453 0.1403 NA Grain Yield
## Mutare_2018.Plant lodging 0.4537 1.2020 0.2740 NA Plant lodging
## Mutare_2018.Plant Height 587.3138 245.3152 0.7054 NA Plant Height
## Mutare_2018.Grain Yield 619789.7743 361342.1945 0.6317 NA Grain Yield
## Mutare_2019.Plant lodging NA NA NA NA Plant lodging
## Mutare_2019.Plant Height NA NA NA NA Plant Height
## Mutare_2019.Grain Yield 208788.9493 204126.3859 0.5056 NA Grain Yield
## Mutare_2020.Plant lodging NA NA NA NA Plant lodging
## Mutare_2020.Plant Height NA NA NA NA Plant Height
## Mutare_2020.Grain Yield 313346.1097 276250.0282 0.5315 NA Grain Yield
## Mutare_2021.Plant lodging 0.0000 3.0276 0.0000 NA Plant lodging
## Mutare_2021.Plant Height NA NA NA NA Plant Height
## Mutare_2021.Grain Yield 99592.7483 618664.7397 0.1387 NA Grain Yield
## Shamva_2018.Plant lodging 0.0782 0.3743 0.1727 NA Plant lodging
## Shamva_2018.Plant Height 493.6128 83.9654 0.8546 NA Plant Height
## Shamva_2018.Grain Yield 43713.4691 150712.6572 0.2248 NA Grain Yield
## Shamva_2019.Plant lodging 0.0954 0.4834 0.1649 NA Plant lodging
## Shamva_2019.Plant Height 174.9860 1440.9164 0.1083 NA Plant Height
## Shamva_2019.Grain Yield 48971.9856 207604.4514 0.1909 NA Grain Yield
## Shamva_2022.Plant lodging 0.0107 0.0218 0.3294 NA Plant lodging
## Shamva_2022.Plant Height 32.7946 114.7649 0.2222 NA Plant Height
## Shamva_2022.Grain Yield 15230.5767 212704.1398 0.0668 NA Grain Yield
## Stapleford_2018.Plant lodging 0.0000 0.2778 0.0000 NA Plant lodging
## Stapleford_2018.Plant Height 353.5130 232.4261 0.6033 NA Plant Height
## Stapleford_2018.Grain Yield 649899.6569 311598.1297 0.6759 NA Grain Yield
## Stapleford_2019.Plant lodging 0.7366 0.3767 0.6616 NA Plant lodging
## Stapleford_2019.Plant Height 81.0351 54.6182 0.5974 NA Plant Height
## Stapleford_2019.Grain Yield 171198.8059 168382.2817 0.5041 NA Grain Yield
## Stapleford_2020.Plant lodging NA NA NA NA Plant lodging
## Stapleford_2020.Plant Height 99.5641 177.8913 0.3588 NA Plant Height
## Stapleford_2020.Grain Yield 653844.4623 442580.0255 0.5963 NA Grain Yield
## Stapleford_2022.Plant lodging 1.1660 0.9009 0.5641 NA Plant lodging
## Stapleford_2022.Plant Height 116.3675 113.0435 0.5072 NA Plant Height
## Stapleford_2022.Grain Yield 311759.2963 301908.6413 0.5080 NA Grain Yield
## Stapleford_2023.Plant lodging NA NA NA NA Plant lodging
## Stapleford_2023.Plant Height 163.7645 137.8405 0.5430 NA Plant Height
## Stapleford_2023.Grain Yield 84997.1803 319570.5776 0.2101 NA Grain Yield
## env geno_sig
## Banket_2018.Plant lodging Banket_2018 *
## Banket_2018.Plant Height Banket_2018 *
## Banket_2018.Grain Yield Banket_2018 *
## Banket_2019.Plant lodging Banket_2019 <NA>
## Banket_2019.Plant Height Banket_2019 <NA>
## Banket_2019.Grain Yield Banket_2019 *
## Banket_2020.Plant lodging Banket_2020 <NA>
## Banket_2020.Plant Height Banket_2020 *
## Banket_2020.Grain Yield Banket_2020 *
## Bindura_2018.Plant lodging Bindura_2018 *
## Bindura_2018.Plant Height Bindura_2018 *
## Bindura_2018.Grain Yield Bindura_2018 *
## Bindura_2020.Plant lodging Bindura_2020 <NA>
## Bindura_2020.Plant Height Bindura_2020 ns
## Bindura_2020.Grain Yield Bindura_2020 *
## Bindura_2021.Plant lodging Bindura_2021 <NA>
## Bindura_2021.Plant Height Bindura_2021 <NA>
## Bindura_2021.Grain Yield Bindura_2021 *
## Chisumbanje_2018.Plant lodging Chisumbanje_2018 *
## Chisumbanje_2018.Plant Height Chisumbanje_2018 *
## Chisumbanje_2018.Grain Yield Chisumbanje_2018 *
## Chisumbanje_2019.Plant lodging Chisumbanje_2019 *
## Chisumbanje_2019.Plant Height Chisumbanje_2019 *
## Chisumbanje_2019.Grain Yield Chisumbanje_2019 ns
## Goromonzi_2018.Plant lodging Goromonzi_2018 <NA>
## Goromonzi_2018.Plant Height Goromonzi_2018 *
## Goromonzi_2018.Grain Yield Goromonzi_2018 *
## Goromonzi_2019.Plant lodging Goromonzi_2019 *
## Goromonzi_2019.Plant Height Goromonzi_2019 *
## Goromonzi_2019.Grain Yield Goromonzi_2019 *
## Goromonzi_2020.Plant lodging Goromonzi_2020 <NA>
## Goromonzi_2020.Plant Height Goromonzi_2020 *
## Goromonzi_2020.Grain Yield Goromonzi_2020 *
## Goromonzi_2021.Plant lodging Goromonzi_2021 <NA>
## Goromonzi_2021.Plant Height Goromonzi_2021 *
## Goromonzi_2021.Grain Yield Goromonzi_2021 *
## Goromonzi_2022.Plant lodging Goromonzi_2022 *
## Goromonzi_2022.Plant Height Goromonzi_2022 *
## Goromonzi_2022.Grain Yield Goromonzi_2022 *
## Goromonzi_2023.Plant lodging Goromonzi_2023 <NA>
## Goromonzi_2023.Plant Height Goromonzi_2023 ns
## Goromonzi_2023.Grain Yield Goromonzi_2023 *
## Gwebi_2020.Plant lodging Gwebi_2020 *
## Gwebi_2020.Plant Height Gwebi_2020 *
## Gwebi_2020.Grain Yield Gwebi_2020 ns
## Gwebi_2021.Plant lodging Gwebi_2021 *
## Gwebi_2021.Plant Height Gwebi_2021 *
## Gwebi_2021.Grain Yield Gwebi_2021 *
## Harare_2018.Plant lodging Harare_2018 *
## Harare_2018.Plant Height Harare_2018 *
## Harare_2018.Grain Yield Harare_2018 *
## Harare_2019.Plant lodging Harare_2019 *
## Harare_2019.Plant Height Harare_2019 *
## Harare_2019.Grain Yield Harare_2019 ns
## Harare_2020.Plant lodging Harare_2020 ns
## Harare_2020.Plant Height Harare_2020 *
## Harare_2020.Grain Yield Harare_2020 *
## Harare_2021.Plant lodging Harare_2021 ns
## Harare_2021.Plant Height Harare_2021 *
## Harare_2021.Grain Yield Harare_2021 *
## Harare_2022.Plant lodging Harare_2022 *
## Harare_2022.Plant Height Harare_2022 ns
## Harare_2022.Grain Yield Harare_2022 ns
## Kadoma_2018.Plant lodging Kadoma_2018 *
## Kadoma_2018.Plant Height Kadoma_2018 ns
## Kadoma_2018.Grain Yield Kadoma_2018 *
## Kadoma_2021.Plant lodging Kadoma_2021 ns
## Kadoma_2021.Plant Height Kadoma_2021 ns
## Kadoma_2021.Grain Yield Kadoma_2021 ns
## Mazowe_2022.Plant lodging Mazowe_2022 ns
## Mazowe_2022.Plant Height Mazowe_2022 *
## Mazowe_2022.Grain Yield Mazowe_2022 ns
## Mthampden_2023.Plant lodging Mthampden_2023 *
## Mthampden_2023.Plant Height Mthampden_2023 *
## Mthampden_2023.Grain Yield Mthampden_2023 ns
## Mutare_2018.Plant lodging Mutare_2018 *
## Mutare_2018.Plant Height Mutare_2018 *
## Mutare_2018.Grain Yield Mutare_2018 *
## Mutare_2019.Plant lodging Mutare_2019 <NA>
## Mutare_2019.Plant Height Mutare_2019 <NA>
## Mutare_2019.Grain Yield Mutare_2019 *
## Mutare_2020.Plant lodging Mutare_2020 <NA>
## Mutare_2020.Plant Height Mutare_2020 <NA>
## Mutare_2020.Grain Yield Mutare_2020 *
## Mutare_2021.Plant lodging Mutare_2021 ns
## Mutare_2021.Plant Height Mutare_2021 <NA>
## Mutare_2021.Grain Yield Mutare_2021 ns
## Shamva_2018.Plant lodging Shamva_2018 ns
## Shamva_2018.Plant Height Shamva_2018 *
## Shamva_2018.Grain Yield Shamva_2018 *
## Shamva_2019.Plant lodging Shamva_2019 *
## Shamva_2019.Plant Height Shamva_2019 ns
## Shamva_2019.Grain Yield Shamva_2019 *
## Shamva_2022.Plant lodging Shamva_2022 *
## Shamva_2022.Plant Height Shamva_2022 ns
## Shamva_2022.Grain Yield Shamva_2022 ns
## Stapleford_2018.Plant lodging Stapleford_2018 ns
## Stapleford_2018.Plant Height Stapleford_2018 *
## Stapleford_2018.Grain Yield Stapleford_2018 *
## Stapleford_2019.Plant lodging Stapleford_2019 *
## Stapleford_2019.Plant Height Stapleford_2019 *
## Stapleford_2019.Grain Yield Stapleford_2019 *
## Stapleford_2020.Plant lodging Stapleford_2020 <NA>
## Stapleford_2020.Plant Height Stapleford_2020 *
## Stapleford_2020.Grain Yield Stapleford_2020 *
## Stapleford_2022.Plant lodging Stapleford_2022 *
## Stapleford_2022.Plant Height Stapleford_2022 *
## Stapleford_2022.Grain Yield Stapleford_2022 *
## Stapleford_2023.Plant lodging Stapleford_2023 <NA>
## Stapleford_2023.Plant Height Stapleford_2023 *
## Stapleford_2023.Grain Yield Stapleford_2023 *Bayesian model
# Filter env which has significance to genotype effect
stmt = readRDS(file = "Saves/stmt.rds")
diagnostic = do.call(rbind, lapply(stmt, function(x){
aa = as.data.frame(x$summary) |> rownames_to_column("trait")
aa
})) |> rownames_to_column("env") |>
mutate_at("env", str_replace, "\\..*", "") |>
dplyr::select(env, trait, geno) |>
separate(geno, into = c("genovar", "sig"), sep = " ")
# diagnostic
# dim(pheno)
# dim(pheno_filt)
# unique(pheno_filt$trial)
# levels(pheno_filt$trial)
# inpheno# install.packages("ProbBreed")
library(ProbBreed)
# head(pheno_filt)
i=traits[3]
# Getting the path of your current open file
# current_path = rstudioapi::getActiveDocumentContext()$path
# setwd(dirname(current_path ))
# print( getwd() )
# head(pheno)
for (i in traits) {
message(paste("Analysis to", i))
# print(paste("Analysis to", i))
# cat("====> Trait:", i, fill = TRUE)
pheno_filt <- droplevels(pheno[
which(pheno$trial %in% diagnostic[which(diagnostic$trait == i &
diagnostic$sig == "*"), "env"]), ])
dim(pheno_filt)
mod = bayes_met(data = pheno_filt,
gen = "geno",
loc = "LOCATION_NAME",
repl = c("rep"),
trait = paste(i),
reg = NULL,
year = "CROP_SEASON",
res.het = TRUE,
iter = 4000, cores = 4, chain = 4)
#save(mod,file = paste0("Saves/",i,".rda"))
saveRDS(mod,file = paste0("Saves/",i,"_year.rds"))
rm(mod)
}
rm(list=ls())
mod_PH=readRDS("Saves/Plant Height.rds")
mod_PL=readRDS("Saves/Plant lodging.rds")
mod_GY=readRDS("Saves/Grain Yield.rds")
# setwd("./Saves/PH")
res_PH = prob_sup(extr = outs_PH,
int = .2,
increase = FALSE,
save.df = FALSE,
verbose = FALSE)
# setwd("../PL")
res_PL = prob_sup(extr = outs_PL,
int = .2,
increase = FALSE,
save.df = FALSE,
verbose = FALSE)
# setwd("../GY")
res_GY = prob_sup(extr = outs_GY,
int = .2,
increase = TRUE,
save.df = FALSE,
verbose = FALSE)
# setwd("../")
save(res_PL,file = "Saves/res_PL_year.rda")
save(res_PH,file = "Saves/res_PH_year.rda")
save(res_GY,file = "Saves/res_GY_year.rda")Variances
outs_PH=readRDS("Saves/outs_PH.rds")
outs_PL=readRDS("Saves/outs_PL.rds")
outs_GY=readRDS("Saves/outs_GY.rds")
outs_PH$variances$trait="PH"
outs_PL$variances$trait="PL"
outs_GY$variances$trait="GY"
varcomp=rbind(outs_GY$variances,outs_PH$variances,outs_PL$variances)
library(tidyverse)
head(varcomp)
varcomp |> ggplot(aes(x=effect,y=var, fill=trait)) + geom_col(position = "dodge") + theme(axis.text.x = element_text(angle=90), legend.position = "none") + facet_wrap(~trait,scales="free")outs_PH=readRDS("Saves/outs_PH_year.rds")
outs_PL=readRDS("Saves/outs_PL_year.rds")
outs_GY=readRDS("Saves/outs_GY_year.rds")
outs_PH$variances$trait="Plant Height"
outs_PL$variances$trait="Plant Lodging"
outs_GY$variances$trait="Grain Yield"
varcomp=rbind(outs_GY$variances,outs_PH$variances,outs_PL$variances)
library(tidyverse)
head(varcomp)
vc=varcomp |> ggplot(aes(x=effect,y=var, fill=trait)) + geom_col(position = "dodge") + theme(axis.text.x = element_text(angle=90), legend.position = "none") +
geom_errorbar(aes(ymin = var-sd,ymax=var+sd),width=0.1)+ facet_wrap(~trait,scales="free")
ggsave(plot=vc,device = "pdf",filename ="Plots/varcomp1.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=vc,device = "png",filename ="Plots/varcomp1.png",dpi = "retina",height = 6,width = 8,bg = "white")Density
library(ProbBreed)
library(tidyverse)
outs_PH=readRDS("Saves/outs_PH_year.rds")
outs_PL=readRDS("Saves/outs_PL_year.rds")
outs_GY=readRDS("Saves/outs_GY_year.rds")
outs_PH$ppcheck
outs_PL$ppcheck
outs_GY$ppcheck
# plot(outs_PH, category = "density")
# plot(outs_PH, category = "histogram")
# plot(outs_PL, category = "density")
# plot(outs_PL, category = "histogram")
# plot(outs_GY, category = "density")
# plot(outs_GY, category = "histogram")
# Density: Empirical vs Sampled phenotypic values
PL=plot(outs_PL)+ggtitle("")+ labs(x="Plant Lodging")
PH=plot(outs_PH)+ggtitle("")
GY=plot(outs_GY)+ggtitle("")
library(ggpubr)
plot.density = ggarrange(GY,PH,PL,labels = c("A","B","C"),nrow = 3, common.legend = T)
ggsave(plot=plot.density,device = "pdf",filename ="Plots/plots_density.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.density,device = "png",filename ="Plots/plots_density.png",dpi = "retina",height = 6,width = 8,bg = "white")
Probsup within locations
rm(list=ls())
library(ProbBreed)
library(tidyverse)
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")
df=read.csv("BPSI/bpsi.csv")
gen.sel = df[which(df$sel=="selected"),"gen"]
GY=plot(res_GY, category = "perfo", level = "within")
PH=plot(res_PH, category = "perfo", level = "within")
PL=plot(res_PL, category = "perfo", level = "within")
unique(GY$data$fac)
cols.label=c("Location", "Year")
names(cols.label)=c(unique(GY$data$fac))
labels <- function(x) {
gsub("G_", "V0", x)
}
GY=GY+ theme(axis.text.y = element_text(size=4,color = ifelse(GY$data$gen %in% gen.sel, "#E41A1C", "black"))) + labs( title = "Grain Yield")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)
PH=PH+ theme(axis.text.y = element_text(size=4, color = ifelse(PH$data$gen %in% gen.sel, "#E41A1C", "black"))) + labs( title = "Plant Height")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)
PL=PL+ theme(axis.text.y = element_text(size=4, color = ifelse(PL$data$gen %in% gen.sel, "#E41A1C", "black"))) + labs( title = "Plant Lodging")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)
library(ggpubr)
plot.within = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol=3, common.legend = T)
ggsave(plot=plot.within,device = "pdf",filename ="Plots/plots_perfo_within.pdf",dpi = "retina",height = 8,width = 10)
ggsave(plot=plot.within,device = "png",filename ="Plots/plots_perfo_within.png",dpi = "retina",height = 8,width = 10,bg = "white")
Probsup across locations
library(ProbBreed)
library(tidyverse)
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")
df=read.csv("BPSI/bpsi.csv")
gen.sel = df[which(df$sel=="selected"),"gen"]
PH=plot(res_PH)
PL=plot(res_PL)
GY=plot(res_GY)
labels <- function(x) {
gsub("G_", "V0", x)
}
library(RColorBrewer)
RColorBrewer::brewer.pal(n=8,name = "Spectral")
display.brewer.all(n=8)
PH=PH + geom_point(aes( x = PH$data$ID,y = PH$data$prob,
fill = ifelse(PH$data$ID %in% gen.sel, "#377EB8","#E41A1C" ))
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +
labs(y = "", x = "", title = "Plant Height")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)
PL=PL + geom_point(aes( x = PL$data$ID,y = PL$data$prob,
fill = ifelse(PL$data$ID %in% gen.sel, "#377EB8","#E41A1C" ))
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +
labs(y = "", x = "", title = "Plant Lodging")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)
GY=GY + geom_point(aes( x = GY$data$ID,y = GY$data$prob,
fill = ifelse(GY$data$ID %in% gen.sel, "#377EB8","#E41A1C" ))
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +
labs(y = "", x = "", title = "Grain Yield")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)
library(ggpubr)
plot.perfo = ggarrange(GY,PH,PL,labels = c("A","B","C"),nrow = 3, common.legend = T)
ggsave(plot=plot.perfo,device = "pdf",filename ="Plots/plots_perfo.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.perfo,device = "png",filename ="Plots/plots_perfo.png",dpi = "retina",height = 6,width = 8)
## Probsup across locations
pairwise
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")
df=read.csv("BPSI/bpsi.csv")
gen.sel = df[which(df$sel=="selected"),"gen"]
GY=plot(res_GY, category = "pair_perfo")
PH=plot(res_PH, category = "pair_perfo")
PL=plot(res_PL, category = "pair_perfo")
library(tidyverse)
target2 = gen.sel
# load("margsarq.rda")
# load("b.rda")
# load("a.rda")
a=GY
b=PH
c=PL
a$data=a$data[which(a$data$x %in% target2 ),]
b$data=b$data[which(b$data$x %in% target2 ),]
c$data=c$data[which(c$data$x %in% target2 ),]
a$data$model ="GY"
b$data$model ="PH"
c$data$model = "PL"
dim(a$data)
highlighted_obs <- c(gen.sel)
GY=a +
theme(
panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
axis.text.y = element_text(color = ifelse( a$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
scale_fill_viridis_c(direction = 1, na.value = "white",
limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Grain Yield",y="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)
# b$data=b$data[which(b$data$x %in% target2 & b$data$y %in% target2),]
PH=b +
theme(
panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
axis.text.y = element_text(color = ifelse( b$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
scale_fill_viridis_c(direction = 1, na.value = "white",
limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Plant Heigth", x="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)
PL=c +
theme( panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
axis.text.y = element_text(color = ifelse( c$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
scale_fill_viridis_c(direction = 1, na.value = "white",
limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Plant Lodging",x="",y="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)
plot.pair = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol = 3, common.legend = T)
plot.pair = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol = 3, common.legend = T)
ggsave(plot=plot.pair,device = "pdf",filename ="Plots/plots_pair.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.pair,device = "png",filename ="Plots/plots_pair.png",dpi = "retina",height = 6,width = 8, bg = "white")
# BPSI - w 2,1,1
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")
source("data/bpsi_fun.R")
models= vector("list",length(3))
models[[1]] = res_GY;
models[[2]] = res_PH
models[[3]] = res_PL
models[[3]]$across$perfo <- models[[3]]$across$perfo[-which(models[[3]]$across$perfo$ID=="G_55"),]
names(models) <- c("GY","PH","PL")
BPSI_soy=BPSI(modlist = models,increase =c(TRUE,FALSE,FALSE),omega = c(2,1,1),int = 0.1,save.df = T,verbose = T )## -> Genotypes selected using BPSI## Warning in dir.create(path = paste0(getwd(), "/BPSI")):
## '/Users/tiagobchagas/Documents/UFV/Zimbabwe_trials/BPSI' already exists## Process completed!## ==> Considering an intensity of 10%, here are the selected candidates:## GY PH PL bpsi sel gen
## 1 0.5 4.0 4.0 8.5 selected G_20
## 2 8.0 6.0 4.0 18.0 selected G_3
## 3 3.0 2.0 17.5 22.5 selected G_38
## 4 7.0 5.0 18.0 30.0 selected G_45
## 5 13.0 3.0 15.0 31.0 selected G_4
## 6 7.0 8.0 24.5 39.5 selected G_13
## 7 29.0 12.0 5.0 46.0 selected G_16
## 8 19.0 18.0 14.0 51.0 selected G_1
## 9 2.0 49.0 1.5 52.5 selected G_32
## 10 19.0 21.5 12.0 52.5 selected G_33
## 11 27.0 7.0 19.0 53.0 not_sel G_14
## 12 15.0 38.0 4.5 57.5 not_sel G_21
## 13 10.0 17.0 33.0 60.0 not_sel G_31
## 14 57.0 1.0 2.5 60.5 not_sel G_54
## 15 2.0 52.0 7.0 61.0 not_sel G_50
## 16 17.0 17.5 28.0 62.5 not_sel G_28
## 17 48.0 3.5 11.0 62.5 not_sel G_64
## 18 5.0 47.0 13.0 65.0 not_sel G_93
## 19 13.0 42.0 13.5 68.5 not_sel G_40
## 20 34.0 24.0 11.5 69.5 not_sel G_49
## 21 23.0 25.0 23.0 71.0 not_sel G_60
## 22 37.0 18.0 17.0 72.0 not_sel G_39
## 23 20.5 39.0 13.0 72.5 not_sel G_68
## 24 36.0 29.0 8.0 73.0 not_sel G_8
## 25 24.0 34.0 18.5 76.5 not_sel G_77
## 26 3.0 61.0 20.0 84.0 not_sel G_53
## 27 18.0 24.0 47.0 89.0 not_sel G_36
## 28 21.0 38.0 30.0 89.0 not_sel G_47
## 29 33.0 20.5 38.0 91.5 not_sel G_52
## 30 23.0 30.0 39.0 92.0 not_sel G_29
## 31 9.0 28.0 55.0 92.0 not_sel G_44
## 32 14.0 76.0 2.0 92.0 not_sel G_9
## 33 32.0 9.5 51.0 92.5 not_sel G_11
## 34 43.0 28.0 21.5 92.5 not_sel G_74
## 35 11.0 58.0 25.0 94.0 not_sel G_42
## 36 44.0 13.0 43.5 100.5 not_sel G_58
## 37 12.0 36.0 53.0 101.0 not_sel G_25
## 38 73.0 10.0 24.0 107.0 not_sel G_67
## 39 11.0 76.0 20.5 107.5 not_sel G_46
## 40 68.0 23.0 18.0 109.0 not_sel G_10
## 41 55.0 4.5 50.0 109.5 not_sel G_2
## 42 49.0 16.0 45.0 110.0 not_sel G_30
## 43 77.0 13.0 21.0 111.0 not_sel G_75
## 44 15.0 53.0 44.0 112.0 not_sel G_92
## 45 54.0 5.0 54.0 113.0 not_sel G_17
## 46 10.0 76.0 29.5 115.5 not_sel G_43
## 47 81.0 5.5 29.0 115.5 not_sel G_5
## 48 31.0 25.0 60.0 116.0 not_sel G_89
## 49 33.5 22.0 63.0 118.5 not_sel G_6
## 50 65.0 15.0 40.0 120.0 not_sel G_88
## 51 8.0 76.0 40.0 124.0 not_sel G_34
## 52 30.5 72.0 22.0 124.5 not_sel G_70
## 53 35.0 76.0 16.0 127.0 not_sel G_63
## 54 51.0 76.0 0.5 127.5 not_sel G_27
## 55 47.0 15.5 66.0 128.5 not_sel G_94
## 56 45.0 60.0 24.0 129.0 not_sel G_35
## 57 50.0 76.0 3.0 129.0 not_sel G_51
## 58 20.0 76.0 34.0 130.0 not_sel G_26
## 59 59.0 22.5 52.0 133.5 not_sel G_61
## 60 42.0 38.0 57.0 137.0 not_sel G_91
## 61 39.0 38.0 64.0 141.0 not_sel G_41
## 62 86.0 16.0 41.0 143.0 not_sel G_19
## 63 31.0 40.0 76.0 147.0 not_sel G_84
## 64 80.0 27.0 42.0 149.0 not_sel G_97
## 65 32.0 27.0 90.0 149.0 not_sel G_98
## 66 63.0 23.0 65.0 151.0 not_sel G_72
## 67 38.5 21.0 93.0 152.5 not_sel G_65
## 68 33.0 59.0 62.0 154.0 not_sel G_81
## 69 56.0 64.0 35.0 155.0 not_sel G_85
## 70 12.5 76.0 67.0 155.5 not_sel G_73
## 71 58.0 64.0 39.5 161.5 not_sel G_82
## 72 60.0 64.0 41.0 165.0 not_sel G_90
## 73 35.5 72.0 58.0 165.5 not_sel G_48
## 74 40.5 33.0 93.0 166.5 not_sel G_57
## 75 43.0 69.0 56.0 168.0 not_sel G_87
## 76 26.0 55.0 88.0 169.0 not_sel G_76
## 77 53.0 28.5 89.0 170.5 not_sel G_80
## 78 42.0 64.0 68.0 174.0 not_sel G_12
## 79 34.5 62.0 84.0 180.5 not_sel G_95
## 80 90.0 76.0 15.5 181.5 not_sel G_24
## 81 37.5 70.0 75.0 182.5 not_sel G_18
## 82 84.0 38.0 61.0 183.0 not_sel G_22
## 83 86.0 63.0 35.5 184.5 not_sel G_71
## 84 86.0 25.5 73.0 184.5 not_sel G_83
## 85 70.0 38.0 77.0 185.0 not_sel G_78
## 86 79.0 18.5 90.0 187.5 not_sel G_56
## 87 74.0 70.0 45.0 189.0 not_sel G_96
## 88 81.0 76.0 36.0 193.0 not_sel G_66
## 89 90.0 22.0 81.0 193.0 not_sel G_7
## 90 72.0 36.0 86.0 194.0 not_sel G_86
## 91 45.0 76.0 74.0 195.0 not_sel G_23
## 92 90.0 38.0 69.0 197.0 not_sel G_59
## 93 45.0 76.0 78.0 199.0 not_sel G_37
## 94 38.0 76.0 85.0 199.0 not_sel G_79
## 95 45.0 64.0 93.0 202.0 not_sel G_15
## 96 90.0 76.0 46.5 212.5 not_sel G_69
## 97 45.0 76.0 93.0 214.0 not_sel G_62BPSI plot - w 2,1,1
Highlighting the ranking of probability of superior performance per trait of the selected families.
Ranking of probability of superior performance per trait - w (2,1,1)
## Warning: Could not calculate the predicate for layer 3; ignored
Selected families based on BPSI rank - w (2,1,1)
##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows| GY | PH | PL | bpsi | sel | gen | |
|---|---|---|---|---|---|---|
| G_20 | 0.5 | 4.0 | 4.0 | 8.5 | selected | G_20 |
| G_3 | 8.0 | 6.0 | 4.0 | 18.0 | selected | G_3 |
| G_38 | 3.0 | 2.0 | 17.5 | 22.5 | selected | G_38 |
| G_45 | 7.0 | 5.0 | 18.0 | 30.0 | selected | G_45 |
| G_4 | 13.0 | 3.0 | 15.0 | 31.0 | selected | G_4 |
| G_13 | 7.0 | 8.0 | 24.5 | 39.5 | selected | G_13 |
| G_16 | 29.0 | 12.0 | 5.0 | 46.0 | selected | G_16 |
| G_1 | 19.0 | 18.0 | 14.0 | 51.0 | selected | G_1 |
| G_32 | 2.0 | 49.0 | 1.5 | 52.5 | selected | G_32 |
| G_33 | 19.0 | 21.5 | 12.0 | 52.5 | selected | G_33 |
| G_14 | 27.0 | 7.0 | 19.0 | 53.0 | not_sel | G_14 |
| G_21 | 15.0 | 38.0 | 4.5 | 57.5 | not_sel | G_21 |
| G_31 | 10.0 | 17.0 | 33.0 | 60.0 | not_sel | G_31 |
| G_54 | 57.0 | 1.0 | 2.5 | 60.5 | not_sel | G_54 |
| G_50 | 2.0 | 52.0 | 7.0 | 61.0 | not_sel | G_50 |
| G_28 | 17.0 | 17.5 | 28.0 | 62.5 | not_sel | G_28 |
| G_64 | 48.0 | 3.5 | 11.0 | 62.5 | not_sel | G_64 |
| G_93 | 5.0 | 47.0 | 13.0 | 65.0 | not_sel | G_93 |
| G_40 | 13.0 | 42.0 | 13.5 | 68.5 | not_sel | G_40 |
| G_49 | 34.0 | 24.0 | 11.5 | 69.5 | not_sel | G_49 |
| G_60 | 23.0 | 25.0 | 23.0 | 71.0 | not_sel | G_60 |
| G_39 | 37.0 | 18.0 | 17.0 | 72.0 | not_sel | G_39 |
| G_68 | 20.5 | 39.0 | 13.0 | 72.5 | not_sel | G_68 |
| G_8 | 36.0 | 29.0 | 8.0 | 73.0 | not_sel | G_8 |
| G_77 | 24.0 | 34.0 | 18.5 | 76.5 | not_sel | G_77 |
| G_53 | 3.0 | 61.0 | 20.0 | 84.0 | not_sel | G_53 |
| G_36 | 18.0 | 24.0 | 47.0 | 89.0 | not_sel | G_36 |
| G_47 | 21.0 | 38.0 | 30.0 | 89.0 | not_sel | G_47 |
| G_52 | 33.0 | 20.5 | 38.0 | 91.5 | not_sel | G_52 |
| G_29 | 23.0 | 30.0 | 39.0 | 92.0 | not_sel | G_29 |
| G_44 | 9.0 | 28.0 | 55.0 | 92.0 | not_sel | G_44 |
| G_9 | 14.0 | 76.0 | 2.0 | 92.0 | not_sel | G_9 |
| G_11 | 32.0 | 9.5 | 51.0 | 92.5 | not_sel | G_11 |
| G_74 | 43.0 | 28.0 | 21.5 | 92.5 | not_sel | G_74 |
| G_42 | 11.0 | 58.0 | 25.0 | 94.0 | not_sel | G_42 |
| G_58 | 44.0 | 13.0 | 43.5 | 100.5 | not_sel | G_58 |
| G_25 | 12.0 | 36.0 | 53.0 | 101.0 | not_sel | G_25 |
| G_67 | 73.0 | 10.0 | 24.0 | 107.0 | not_sel | G_67 |
| G_46 | 11.0 | 76.0 | 20.5 | 107.5 | not_sel | G_46 |
| G_10 | 68.0 | 23.0 | 18.0 | 109.0 | not_sel | G_10 |
| G_2 | 55.0 | 4.5 | 50.0 | 109.5 | not_sel | G_2 |
| G_30 | 49.0 | 16.0 | 45.0 | 110.0 | not_sel | G_30 |
| G_75 | 77.0 | 13.0 | 21.0 | 111.0 | not_sel | G_75 |
| G_92 | 15.0 | 53.0 | 44.0 | 112.0 | not_sel | G_92 |
| G_17 | 54.0 | 5.0 | 54.0 | 113.0 | not_sel | G_17 |
| G_43 | 10.0 | 76.0 | 29.5 | 115.5 | not_sel | G_43 |
| G_5 | 81.0 | 5.5 | 29.0 | 115.5 | not_sel | G_5 |
| G_89 | 31.0 | 25.0 | 60.0 | 116.0 | not_sel | G_89 |
| G_6 | 33.5 | 22.0 | 63.0 | 118.5 | not_sel | G_6 |
| G_88 | 65.0 | 15.0 | 40.0 | 120.0 | not_sel | G_88 |
| G_34 | 8.0 | 76.0 | 40.0 | 124.0 | not_sel | G_34 |
| G_70 | 30.5 | 72.0 | 22.0 | 124.5 | not_sel | G_70 |
| G_63 | 35.0 | 76.0 | 16.0 | 127.0 | not_sel | G_63 |
| G_27 | 51.0 | 76.0 | 0.5 | 127.5 | not_sel | G_27 |
| G_94 | 47.0 | 15.5 | 66.0 | 128.5 | not_sel | G_94 |
| G_35 | 45.0 | 60.0 | 24.0 | 129.0 | not_sel | G_35 |
| G_51 | 50.0 | 76.0 | 3.0 | 129.0 | not_sel | G_51 |
| G_26 | 20.0 | 76.0 | 34.0 | 130.0 | not_sel | G_26 |
| G_61 | 59.0 | 22.5 | 52.0 | 133.5 | not_sel | G_61 |
| G_91 | 42.0 | 38.0 | 57.0 | 137.0 | not_sel | G_91 |
| G_41 | 39.0 | 38.0 | 64.0 | 141.0 | not_sel | G_41 |
| G_19 | 86.0 | 16.0 | 41.0 | 143.0 | not_sel | G_19 |
| G_84 | 31.0 | 40.0 | 76.0 | 147.0 | not_sel | G_84 |
| G_97 | 80.0 | 27.0 | 42.0 | 149.0 | not_sel | G_97 |
| G_98 | 32.0 | 27.0 | 90.0 | 149.0 | not_sel | G_98 |
| G_72 | 63.0 | 23.0 | 65.0 | 151.0 | not_sel | G_72 |
| G_65 | 38.5 | 21.0 | 93.0 | 152.5 | not_sel | G_65 |
| G_81 | 33.0 | 59.0 | 62.0 | 154.0 | not_sel | G_81 |
| G_85 | 56.0 | 64.0 | 35.0 | 155.0 | not_sel | G_85 |
| G_73 | 12.5 | 76.0 | 67.0 | 155.5 | not_sel | G_73 |
| G_82 | 58.0 | 64.0 | 39.5 | 161.5 | not_sel | G_82 |
| G_90 | 60.0 | 64.0 | 41.0 | 165.0 | not_sel | G_90 |
| G_48 | 35.5 | 72.0 | 58.0 | 165.5 | not_sel | G_48 |
| G_57 | 40.5 | 33.0 | 93.0 | 166.5 | not_sel | G_57 |
| G_87 | 43.0 | 69.0 | 56.0 | 168.0 | not_sel | G_87 |
| G_76 | 26.0 | 55.0 | 88.0 | 169.0 | not_sel | G_76 |
| G_80 | 53.0 | 28.5 | 89.0 | 170.5 | not_sel | G_80 |
| G_12 | 42.0 | 64.0 | 68.0 | 174.0 | not_sel | G_12 |
| G_95 | 34.5 | 62.0 | 84.0 | 180.5 | not_sel | G_95 |
| G_24 | 90.0 | 76.0 | 15.5 | 181.5 | not_sel | G_24 |
| G_18 | 37.5 | 70.0 | 75.0 | 182.5 | not_sel | G_18 |
| G_22 | 84.0 | 38.0 | 61.0 | 183.0 | not_sel | G_22 |
| G_71 | 86.0 | 63.0 | 35.5 | 184.5 | not_sel | G_71 |
| G_83 | 86.0 | 25.5 | 73.0 | 184.5 | not_sel | G_83 |
| G_78 | 70.0 | 38.0 | 77.0 | 185.0 | not_sel | G_78 |
| G_56 | 79.0 | 18.5 | 90.0 | 187.5 | not_sel | G_56 |
| G_96 | 74.0 | 70.0 | 45.0 | 189.0 | not_sel | G_96 |
| G_66 | 81.0 | 76.0 | 36.0 | 193.0 | not_sel | G_66 |
| G_7 | 90.0 | 22.0 | 81.0 | 193.0 | not_sel | G_7 |
| G_86 | 72.0 | 36.0 | 86.0 | 194.0 | not_sel | G_86 |
| G_23 | 45.0 | 76.0 | 74.0 | 195.0 | not_sel | G_23 |
| G_59 | 90.0 | 38.0 | 69.0 | 197.0 | not_sel | G_59 |
| G_37 | 45.0 | 76.0 | 78.0 | 199.0 | not_sel | G_37 |
| G_79 | 38.0 | 76.0 | 85.0 | 199.0 | not_sel | G_79 |
| G_15 | 45.0 | 64.0 | 93.0 | 202.0 | not_sel | G_15 |
| G_69 | 90.0 | 76.0 | 46.5 | 212.5 | not_sel | G_69 |
| G_62 | 45.0 | 76.0 | 93.0 | 214.0 | not_sel | G_62 |
BPSI - w 3,1,1
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")
source("data/bpsi_fun.R")
models= vector("list",length(3))
models[[1]] = res_GY;
models[[2]] = res_PH
models[[3]] = res_PL
models[[3]]$across$perfo <- models[[3]]$across$perfo[-which(models[[3]]$across$perfo$ID=="G_55"),]
names(models) <- c("GY","PH","PL")
BPSI_soy=BPSI(modlist = models,increase =c(TRUE,FALSE,FALSE),omega = c(3,1,1),int = 0.1,save.df = T,verbose = T )## -> Genotypes selected using BPSI## Warning in dir.create(path = paste0(getwd(), "/BPSI")):
## '/Users/tiagobchagas/Documents/UFV/Zimbabwe_trials/BPSI' already exists## Process completed!## ==> Considering an intensity of 10%, here are the selected candidates:## GY PH PL bpsi sel gen
## 1 0.3333333 4.000000 4.0000000 8.333333 selected G_20
## 2 3.0000000 2.000000 11.6666667 16.666667 selected G_38
## 3 8.0000000 6.000000 2.6666667 16.666667 selected G_3
## 4 8.6666667 3.000000 15.0000000 26.666667 selected G_4
## 5 4.6666667 5.000000 18.0000000 27.666667 selected G_45
## 6 7.0000000 8.000000 16.3333333 31.333333 selected G_13
## 7 29.0000000 12.000000 3.3333333 44.333333 selected G_16
## 8 12.6666667 18.000000 14.0000000 44.666667 selected G_1
## 9 19.0000000 14.333333 12.0000000 45.333333 selected G_33
## 10 27.0000000 4.666667 19.0000000 50.666667 selected G_14
## 11 2.0000000 49.000000 1.0000000 52.000000 not_sel G_32
## 12 15.0000000 38.000000 3.0000000 56.000000 not_sel G_21
## 13 17.0000000 11.666667 28.0000000 56.666667 not_sel G_28
## 14 6.6666667 17.000000 33.0000000 56.666667 not_sel G_31
## 15 57.0000000 1.000000 1.6666667 59.666667 not_sel G_54
## 16 1.3333333 52.000000 7.0000000 60.333333 not_sel G_50
## 17 5.0000000 47.000000 8.6666667 60.666667 not_sel G_93
## 18 48.0000000 2.333333 11.0000000 61.333333 not_sel G_64
## 19 23.0000000 25.000000 15.3333333 63.333333 not_sel G_60
## 20 13.0000000 42.000000 9.0000000 64.000000 not_sel G_40
## 21 13.6666667 39.000000 13.0000000 65.666667 not_sel G_68
## 22 34.0000000 24.000000 7.6666667 65.666667 not_sel G_49
## 23 37.0000000 12.000000 17.0000000 66.000000 not_sel G_39
## 24 24.0000000 34.000000 12.3333333 70.333333 not_sel G_77
## 25 36.0000000 29.000000 5.3333333 70.333333 not_sel G_8
## 26 21.0000000 25.333333 30.0000000 76.333333 not_sel G_47
## 27 18.0000000 16.000000 47.0000000 81.000000 not_sel G_36
## 28 9.0000000 18.666667 55.0000000 82.666667 not_sel G_44
## 29 2.0000000 61.000000 20.0000000 83.000000 not_sel G_53
## 30 15.3333333 30.000000 39.0000000 84.333333 not_sel G_29
## 31 33.0000000 13.666667 38.0000000 84.666667 not_sel G_52
## 32 43.0000000 28.000000 14.3333333 85.333333 not_sel G_74
## 33 44.0000000 13.000000 29.0000000 86.000000 not_sel G_58
## 34 9.3333333 76.000000 2.0000000 87.333333 not_sel G_9
## 35 12.0000000 24.000000 53.0000000 89.000000 not_sel G_25
## 36 32.0000000 6.333333 51.0000000 89.333333 not_sel G_11
## 37 7.3333333 58.000000 25.0000000 90.333333 not_sel G_42
## 38 11.0000000 76.000000 13.6666667 100.666667 not_sel G_46
## 39 68.0000000 23.000000 12.0000000 103.000000 not_sel G_10
## 40 73.0000000 6.666667 24.0000000 103.666667 not_sel G_67
## 41 49.0000000 10.666667 45.0000000 104.666667 not_sel G_30
## 42 10.0000000 76.000000 19.6666667 105.666667 not_sel G_43
## 43 77.0000000 8.666667 21.0000000 106.666667 not_sel G_75
## 44 65.0000000 15.000000 26.6666667 106.666667 not_sel G_88
## 45 10.0000000 53.000000 44.0000000 107.000000 not_sel G_92
## 46 22.3333333 22.000000 63.0000000 107.333333 not_sel G_6
## 47 31.0000000 16.666667 60.0000000 107.666667 not_sel G_89
## 48 55.0000000 3.000000 50.0000000 108.000000 not_sel G_2
## 49 54.0000000 3.333333 54.0000000 111.333333 not_sel G_17
## 50 81.0000000 3.666667 29.0000000 113.666667 not_sel G_5
## 51 20.3333333 72.000000 22.0000000 114.333333 not_sel G_70
## 52 45.0000000 60.000000 16.0000000 121.000000 not_sel G_35
## 53 5.3333333 76.000000 40.0000000 121.333333 not_sel G_34
## 54 35.0000000 76.000000 10.6666667 121.666667 not_sel G_63
## 55 13.3333333 76.000000 34.0000000 123.333333 not_sel G_26
## 56 47.0000000 10.333333 66.0000000 123.333333 not_sel G_94
## 57 42.0000000 25.333333 57.0000000 124.333333 not_sel G_91
## 58 59.0000000 15.000000 52.0000000 126.000000 not_sel G_61
## 59 51.0000000 76.000000 0.3333333 127.333333 not_sel G_27
## 60 50.0000000 76.000000 2.0000000 128.000000 not_sel G_51
## 61 39.0000000 25.333333 64.0000000 128.333333 not_sel G_41
## 62 86.0000000 16.000000 27.3333333 129.333333 not_sel G_19
## 63 20.6666667 40.000000 76.0000000 136.666667 not_sel G_84
## 64 21.3333333 27.000000 90.0000000 138.333333 not_sel G_98
## 65 25.6666667 21.000000 93.0000000 139.666667 not_sel G_65
## 66 80.0000000 18.000000 42.0000000 140.000000 not_sel G_97
## 67 22.0000000 59.000000 62.0000000 143.000000 not_sel G_81
## 68 63.0000000 15.333333 65.0000000 143.333333 not_sel G_72
## 69 56.0000000 64.000000 23.3333333 143.333333 not_sel G_85
## 70 58.0000000 64.000000 26.3333333 148.333333 not_sel G_82
## 71 8.3333333 76.000000 67.0000000 151.333333 not_sel G_73
## 72 60.0000000 64.000000 27.3333333 151.333333 not_sel G_90
## 73 27.0000000 33.000000 93.0000000 153.000000 not_sel G_57
## 74 23.6666667 72.000000 58.0000000 153.666667 not_sel G_48
## 75 28.6666667 69.000000 56.0000000 153.666667 not_sel G_87
## 76 28.0000000 64.000000 68.0000000 160.000000 not_sel G_12
## 77 17.3333333 55.000000 88.0000000 160.333333 not_sel G_76
## 78 53.0000000 19.000000 89.0000000 161.000000 not_sel G_80
## 79 23.0000000 62.000000 84.0000000 169.000000 not_sel G_95
## 80 25.0000000 70.000000 75.0000000 170.000000 not_sel G_18
## 81 84.0000000 25.333333 61.0000000 170.333333 not_sel G_22
## 82 70.0000000 25.333333 77.0000000 172.333333 not_sel G_78
## 83 86.0000000 63.000000 23.6666667 172.666667 not_sel G_71
## 84 74.0000000 70.000000 30.0000000 174.000000 not_sel G_96
## 85 86.0000000 17.000000 73.0000000 176.000000 not_sel G_83
## 86 90.0000000 76.000000 10.3333333 176.333333 not_sel G_24
## 87 30.0000000 76.000000 74.0000000 180.000000 not_sel G_23
## 88 81.0000000 76.000000 24.0000000 181.000000 not_sel G_66
## 89 79.0000000 12.333333 90.0000000 181.333333 not_sel G_56
## 90 72.0000000 24.000000 86.0000000 182.000000 not_sel G_86
## 91 30.0000000 76.000000 78.0000000 184.000000 not_sel G_37
## 92 90.0000000 25.333333 69.0000000 184.333333 not_sel G_59
## 93 90.0000000 14.666667 81.0000000 185.666667 not_sel G_7
## 94 25.3333333 76.000000 85.0000000 186.333333 not_sel G_79
## 95 30.0000000 64.000000 93.0000000 187.000000 not_sel G_15
## 96 90.0000000 76.000000 31.0000000 197.000000 not_sel G_69
## 97 30.0000000 76.000000 93.0000000 199.000000 not_sel G_62BPSI plot - w 3,1,1
Highlighting the ranking of probability of superior performance per trait of the selected families.
Ranking of probability of superior performance per trait - w (3,1,1)
## Warning: Could not calculate the predicate for layer 3; ignored
Selected families based on BPSI rank - w (3,1,1)
| GY | PH | PL | bpsi | sel | gen | |
|---|---|---|---|---|---|---|
| G_20 | 0.3333333 | 4.000000 | 4.0000000 | 8.333333 | selected | G_20 |
| G_38 | 3.0000000 | 2.000000 | 11.6666667 | 16.666667 | selected | G_38 |
| G_3 | 8.0000000 | 6.000000 | 2.6666667 | 16.666667 | selected | G_3 |
| G_4 | 8.6666667 | 3.000000 | 15.0000000 | 26.666667 | selected | G_4 |
| G_45 | 4.6666667 | 5.000000 | 18.0000000 | 27.666667 | selected | G_45 |
| G_13 | 7.0000000 | 8.000000 | 16.3333333 | 31.333333 | selected | G_13 |
| G_16 | 29.0000000 | 12.000000 | 3.3333333 | 44.333333 | selected | G_16 |
| G_1 | 12.6666667 | 18.000000 | 14.0000000 | 44.666667 | selected | G_1 |
| G_33 | 19.0000000 | 14.333333 | 12.0000000 | 45.333333 | selected | G_33 |
| G_14 | 27.0000000 | 4.666667 | 19.0000000 | 50.666667 | selected | G_14 |
| G_32 | 2.0000000 | 49.000000 | 1.0000000 | 52.000000 | not_sel | G_32 |
| G_21 | 15.0000000 | 38.000000 | 3.0000000 | 56.000000 | not_sel | G_21 |
| G_28 | 17.0000000 | 11.666667 | 28.0000000 | 56.666667 | not_sel | G_28 |
| G_31 | 6.6666667 | 17.000000 | 33.0000000 | 56.666667 | not_sel | G_31 |
| G_54 | 57.0000000 | 1.000000 | 1.6666667 | 59.666667 | not_sel | G_54 |
| G_50 | 1.3333333 | 52.000000 | 7.0000000 | 60.333333 | not_sel | G_50 |
| G_93 | 5.0000000 | 47.000000 | 8.6666667 | 60.666667 | not_sel | G_93 |
| G_64 | 48.0000000 | 2.333333 | 11.0000000 | 61.333333 | not_sel | G_64 |
| G_60 | 23.0000000 | 25.000000 | 15.3333333 | 63.333333 | not_sel | G_60 |
| G_40 | 13.0000000 | 42.000000 | 9.0000000 | 64.000000 | not_sel | G_40 |
| G_68 | 13.6666667 | 39.000000 | 13.0000000 | 65.666667 | not_sel | G_68 |
| G_49 | 34.0000000 | 24.000000 | 7.6666667 | 65.666667 | not_sel | G_49 |
| G_39 | 37.0000000 | 12.000000 | 17.0000000 | 66.000000 | not_sel | G_39 |
| G_77 | 24.0000000 | 34.000000 | 12.3333333 | 70.333333 | not_sel | G_77 |
| G_8 | 36.0000000 | 29.000000 | 5.3333333 | 70.333333 | not_sel | G_8 |
| G_47 | 21.0000000 | 25.333333 | 30.0000000 | 76.333333 | not_sel | G_47 |
| G_36 | 18.0000000 | 16.000000 | 47.0000000 | 81.000000 | not_sel | G_36 |
| G_44 | 9.0000000 | 18.666667 | 55.0000000 | 82.666667 | not_sel | G_44 |
| G_53 | 2.0000000 | 61.000000 | 20.0000000 | 83.000000 | not_sel | G_53 |
| G_29 | 15.3333333 | 30.000000 | 39.0000000 | 84.333333 | not_sel | G_29 |
| G_52 | 33.0000000 | 13.666667 | 38.0000000 | 84.666667 | not_sel | G_52 |
| G_74 | 43.0000000 | 28.000000 | 14.3333333 | 85.333333 | not_sel | G_74 |
| G_58 | 44.0000000 | 13.000000 | 29.0000000 | 86.000000 | not_sel | G_58 |
| G_9 | 9.3333333 | 76.000000 | 2.0000000 | 87.333333 | not_sel | G_9 |
| G_25 | 12.0000000 | 24.000000 | 53.0000000 | 89.000000 | not_sel | G_25 |
| G_11 | 32.0000000 | 6.333333 | 51.0000000 | 89.333333 | not_sel | G_11 |
| G_42 | 7.3333333 | 58.000000 | 25.0000000 | 90.333333 | not_sel | G_42 |
| G_46 | 11.0000000 | 76.000000 | 13.6666667 | 100.666667 | not_sel | G_46 |
| G_10 | 68.0000000 | 23.000000 | 12.0000000 | 103.000000 | not_sel | G_10 |
| G_67 | 73.0000000 | 6.666667 | 24.0000000 | 103.666667 | not_sel | G_67 |
| G_30 | 49.0000000 | 10.666667 | 45.0000000 | 104.666667 | not_sel | G_30 |
| G_43 | 10.0000000 | 76.000000 | 19.6666667 | 105.666667 | not_sel | G_43 |
| G_75 | 77.0000000 | 8.666667 | 21.0000000 | 106.666667 | not_sel | G_75 |
| G_88 | 65.0000000 | 15.000000 | 26.6666667 | 106.666667 | not_sel | G_88 |
| G_92 | 10.0000000 | 53.000000 | 44.0000000 | 107.000000 | not_sel | G_92 |
| G_6 | 22.3333333 | 22.000000 | 63.0000000 | 107.333333 | not_sel | G_6 |
| G_89 | 31.0000000 | 16.666667 | 60.0000000 | 107.666667 | not_sel | G_89 |
| G_2 | 55.0000000 | 3.000000 | 50.0000000 | 108.000000 | not_sel | G_2 |
| G_17 | 54.0000000 | 3.333333 | 54.0000000 | 111.333333 | not_sel | G_17 |
| G_5 | 81.0000000 | 3.666667 | 29.0000000 | 113.666667 | not_sel | G_5 |
| G_70 | 20.3333333 | 72.000000 | 22.0000000 | 114.333333 | not_sel | G_70 |
| G_35 | 45.0000000 | 60.000000 | 16.0000000 | 121.000000 | not_sel | G_35 |
| G_34 | 5.3333333 | 76.000000 | 40.0000000 | 121.333333 | not_sel | G_34 |
| G_63 | 35.0000000 | 76.000000 | 10.6666667 | 121.666667 | not_sel | G_63 |
| G_26 | 13.3333333 | 76.000000 | 34.0000000 | 123.333333 | not_sel | G_26 |
| G_94 | 47.0000000 | 10.333333 | 66.0000000 | 123.333333 | not_sel | G_94 |
| G_91 | 42.0000000 | 25.333333 | 57.0000000 | 124.333333 | not_sel | G_91 |
| G_61 | 59.0000000 | 15.000000 | 52.0000000 | 126.000000 | not_sel | G_61 |
| G_27 | 51.0000000 | 76.000000 | 0.3333333 | 127.333333 | not_sel | G_27 |
| G_51 | 50.0000000 | 76.000000 | 2.0000000 | 128.000000 | not_sel | G_51 |
| G_41 | 39.0000000 | 25.333333 | 64.0000000 | 128.333333 | not_sel | G_41 |
| G_19 | 86.0000000 | 16.000000 | 27.3333333 | 129.333333 | not_sel | G_19 |
| G_84 | 20.6666667 | 40.000000 | 76.0000000 | 136.666667 | not_sel | G_84 |
| G_98 | 21.3333333 | 27.000000 | 90.0000000 | 138.333333 | not_sel | G_98 |
| G_65 | 25.6666667 | 21.000000 | 93.0000000 | 139.666667 | not_sel | G_65 |
| G_97 | 80.0000000 | 18.000000 | 42.0000000 | 140.000000 | not_sel | G_97 |
| G_81 | 22.0000000 | 59.000000 | 62.0000000 | 143.000000 | not_sel | G_81 |
| G_72 | 63.0000000 | 15.333333 | 65.0000000 | 143.333333 | not_sel | G_72 |
| G_85 | 56.0000000 | 64.000000 | 23.3333333 | 143.333333 | not_sel | G_85 |
| G_82 | 58.0000000 | 64.000000 | 26.3333333 | 148.333333 | not_sel | G_82 |
| G_73 | 8.3333333 | 76.000000 | 67.0000000 | 151.333333 | not_sel | G_73 |
| G_90 | 60.0000000 | 64.000000 | 27.3333333 | 151.333333 | not_sel | G_90 |
| G_57 | 27.0000000 | 33.000000 | 93.0000000 | 153.000000 | not_sel | G_57 |
| G_48 | 23.6666667 | 72.000000 | 58.0000000 | 153.666667 | not_sel | G_48 |
| G_87 | 28.6666667 | 69.000000 | 56.0000000 | 153.666667 | not_sel | G_87 |
| G_12 | 28.0000000 | 64.000000 | 68.0000000 | 160.000000 | not_sel | G_12 |
| G_76 | 17.3333333 | 55.000000 | 88.0000000 | 160.333333 | not_sel | G_76 |
| G_80 | 53.0000000 | 19.000000 | 89.0000000 | 161.000000 | not_sel | G_80 |
| G_95 | 23.0000000 | 62.000000 | 84.0000000 | 169.000000 | not_sel | G_95 |
| G_18 | 25.0000000 | 70.000000 | 75.0000000 | 170.000000 | not_sel | G_18 |
| G_22 | 84.0000000 | 25.333333 | 61.0000000 | 170.333333 | not_sel | G_22 |
| G_78 | 70.0000000 | 25.333333 | 77.0000000 | 172.333333 | not_sel | G_78 |
| G_71 | 86.0000000 | 63.000000 | 23.6666667 | 172.666667 | not_sel | G_71 |
| G_96 | 74.0000000 | 70.000000 | 30.0000000 | 174.000000 | not_sel | G_96 |
| G_83 | 86.0000000 | 17.000000 | 73.0000000 | 176.000000 | not_sel | G_83 |
| G_24 | 90.0000000 | 76.000000 | 10.3333333 | 176.333333 | not_sel | G_24 |
| G_23 | 30.0000000 | 76.000000 | 74.0000000 | 180.000000 | not_sel | G_23 |
| G_66 | 81.0000000 | 76.000000 | 24.0000000 | 181.000000 | not_sel | G_66 |
| G_56 | 79.0000000 | 12.333333 | 90.0000000 | 181.333333 | not_sel | G_56 |
| G_86 | 72.0000000 | 24.000000 | 86.0000000 | 182.000000 | not_sel | G_86 |
| G_37 | 30.0000000 | 76.000000 | 78.0000000 | 184.000000 | not_sel | G_37 |
| G_59 | 90.0000000 | 25.333333 | 69.0000000 | 184.333333 | not_sel | G_59 |
| G_7 | 90.0000000 | 14.666667 | 81.0000000 | 185.666667 | not_sel | G_7 |
| G_79 | 25.3333333 | 76.000000 | 85.0000000 | 186.333333 | not_sel | G_79 |
| G_15 | 30.0000000 | 64.000000 | 93.0000000 | 187.000000 | not_sel | G_15 |
| G_69 | 90.0000000 | 76.000000 | 31.0000000 | 197.000000 | not_sel | G_69 |
| G_62 | 30.0000000 | 76.000000 | 93.0000000 | 199.000000 | not_sel | G_62 |